Effects of gap anisotropy upon the electronic structure around a superconducting 

vortex 
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An isolated single vortex is considered within the frame- 
work of the quasiclassical theory. The local density of states 
around a vortex is calculated in a clean type II superconduc- 
tor with an anisotropy. The anisotropy of a superconducting 
energy gap is crucial for bound states around a vortex. A 
characteristic structure of the local density of states, observed 
in the layered hexagonal superconductor 2_ff-NbSe2 by scan- 
ning tunneling microscopy (STM), is well reproduced if one 
assumes an anisotropic s-wave gap in the hexagonal plane. 
The local density of states (or the bound states) around the 
vortex is interpreted in terms of quasiparticle trajectories to 
facilitate an understanding of the rich electronic structure ob- 
served in STM experiments. It is pointed out that further fine 
structures and extra peaks in the local density of states should 
be observed by STM. 

PACS number(s): 74.60.Ec, 61.16.Ch, 74.25.Jb 



I. INTRODUCTION 

Recently, the existence of an anisotropy of a supercon- 
ducting energy gap has attracted a great deal of atten- 
tion in various superconductors such as heavy Fermion, 
organic, and high T c compounds. On the other hand, 
the electronic structure around vortices is a fundamental 
problem on the physics of both conventional and uncon- 
ventional superconductors. In this paper, we discuss ef- 
fects of the gap anisotropy upon the electronic structure 
around a vortex, i.e., the bound states around an isolated 
vortex in clean type II superconductors. 

Theoretically, the bound states around a vortex was. 
discussed in 1964 by Caroli, de Gennes, and Matricon,EJ 
who considered a single vortex in an isotropic s-wave su- 
perconductor. After this work, several theorists stud- 
ied the electronic structure around vortices.HLJ Exper- 
imentally, however, until the following success by Hess 
et al, there had existed for a long time no experi- 
ments which could directly study the electronic structure 
around vortices.B 

In 1989, a novel experimental method, scanning tun- 
neling microscopy (STM), opened up a way to the study 
of the electronicuSjtructure around vortices in type II 
supjerconductorsBQ Using the STM method, Hess et 
al.u succeeded in measuring spatially resolved excitation 
spectra, i.e., local density of states (LDOS) around a 
vortex. They investigated the bound states around a 
vortex in the layered hexagonal compound 2i7-NbSe2 



(T c =7.2 K), and found a striking zero-bias peak at the 
vortex center. The same peak and its collapse upon sub- 
stituting Ta for Nb as impurities in NbSe2 were also ob- 
served by Renner et alu Several new theoreticaLstud- 
ies of the electronic structure around a vortexE2rb3 were 
promoted by the success of the-STMuexperiment by Hess 
et alu Some of these theoriesE3Ejlill predicted that the 
zero-bias peak should split into two, i.e., into positive and 
negative bias voltage peaks, if spectra are taken at some 
distance from the vortex center (see, for instance, Fig. 3 
in Ref. |l^) . This splitting indicates that quasiparticles of 
the vortex bound states with finite angular momentum 
are distributed circularly, and circulate farther away from 
the core center as they have higher energy. The predicted 
splitting was actually confirmed in an experiment O 

However, a mystery also emerged. In the above exper- 
iment, Hess et aLlIa not only confirmed the splitting, but 
also found that the LDOS around the vortex was shaped 
like "star" at a fixed energy and its orientation was de- 
pendent on the energy, i.e., the sixfold star shape rotates 
as the bias voltage varies. (Fig. 4 in Ref. |9|)-Soon after 
this observation was made, Gygi and SchliiterE] proposed 
an explanation for this 30° rotation of the star-shaped 
LDOS. On the basis of a sixfold perturbation, they ex- 
plained that the two states, i.e., the lower and higher 
energy stacs were interpreted as bonding or antibond- 
ing states.Ej Although they explained certain aspects of 
the observation, the following features_of-thfi-star-shaped 
LDOS observed in later experimentscJ&QEj could not 
be sufficiently understood by this perturbation scheme. 
In the intermediate energy, a r^ffiA' of the star splits into 
a pair of nearly parallel rays.EilE3 (Fig. 1 in Ref. |2l] or 
Fig. 1 in Ref. g2|.) The zero-bias peak in the spectral 
evolution along a radial line from the vortex center does 
not split into two subpeaks observedj-ip-|the earlier exper- 
iment, but into three or more ones.LrnJ (Figs. 9 and 10 
in Ref. 0, or Fig. 6 in Ref. ^|.) 

Specifically, the jd^.a£terjstic features of the LDOS 
observed in NbSe2c3tijO'l3o are summarized in detail 
as follows, when the magnetic field H is applied perpen- 
dicular to the hexagonal plane: (1) The LDOS for quasi- 
particle excitations-,has a sixfold star shape centered at 
the vortex centerJia (2) The orientation of this star de- 
pends on the energy. At zero bias, a ray of the star 
extends away from the a axis in the hexagonal plane of 
NbSe2- Upon increasing the bias voltage, the star ro- 
tates by 30°. (3) In the intermediate bias voltage, a ray 
splits intO|-a.pair of nearly parallel rays, keeping its direc- 
tion fixedEJ't3 (4) In the spectral evolution which crosses 
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the vortex center, there exist inner peaks in addition to 
the outer peaks which evolve from the zero bias peak at 
the vortex centej.-into the bulk BCS like gap edges far 
from the vortex. 0E3 The inner peaks vary with the angle 
of the direction in which the spectral evolution is taken. 
These important and interesting observations (l)-(4) re- 
main unexplained. 

Quite recently, motivated by a possibility of a d-wave 
supeixtmductivity in high T c cuprates, Schopohl and 
MakiEJ'Ej studied the electronic structure around a vor- 
tex in a d-wave superconductor. On pthe— basis of the 
quasiclassical Green's function theory,c3c3 the LDOS 
around a single vortex was calculated in a superconduc- 
tor with a d-wave energy gap. They found that the LDOS 
exhibits a characteristic fourfold structure in the d-wave 
gap case, which is contrasted with the isotropic s-wave 
gap case (a circularly symmetric LD0S)O'E£l A gradual 
45° rotation of this fourfold LDOS as the [energy changes 
was later reported by the present authors. Ej We note. that 
this rotation is similar to that observed in NbSe2 £9 

In the context described above, we have investigated 
the electronic structure around the vortex observed in 
NbSe 2 . We find that the rich sttnctuxa p£ the LDOS ob- 
served in the STM experimentaljOL3lJa results mainly 
from a superconducting gap anisotropy. Assuming an 
anisotropic s-wave gap analogous to the d-wave one, we 
are able to obtain results favorably comparable with the 
experiments. In a previous short paperj^d we enumer- 
ated the following items as the possible origin of the rich 
structure of the LDOS: (a) the effect of an anisotropic 
superconducting energy gap, (b) the effect of nearest- 
neighbor vortices, i.e., the effect of the vortex lattice, 
and (c) the effect of the anisotropic density of states at 
the Fermi surface. It is the purpose of the present paper 
to discuss the gap effect (a) in more detail. As for the 
item (b), i.e., the vortex lattice effect, we gave a detailed 
report in Ref. [32]. A detailed report for (c) the effect of 
the anisotropic density of states at the Fermi surface is 
given elsewhere £3 

To date, NbSe2 has been the only compound in which 
the electronic structure around vortices was thoroughly 
investigated by STM. In this paper, we concentrate our 
attention on the LDOS observed in NbSe2 as a typi- 
cal example of a type II superconductor. However, the 
essence of the present considerations is equally applicable 
to other type II superconductors in general. The LDOS 
around a vortex reflects the internal electronic structure 
of the vortex, and an understanding of this structure is 
important in elucidating dynamical properties of vortices 
as well as static ones. 

We consider the case of an isolated static vortex un- 
der a magnetic field applied parallel to the c axis (or z 
axis). We restrict ourselves to a two dimensional system, 
i.e., assume a two dimensional Fermi surface neglecting 
a small warping of the Fermi surface along the c axis, 
which ^appropriate to layered superconductors such as 
NbSe 2 MS 

In Sec. P, we describe the quasiclassical theory we use 



for the study of the vortex. Section III is devoted to 



the calculations of the LDOS around a vortex under the 



influence of the gap anisotropy. In Sec. IV, we interpret 



the resultant LDOS in terms of quasiparticle trajectories. 
The summary and discussions are given in Sec. 0. 



II. THE QUASICLASSICAL THEORY 

To investigate the LDOS around a vortex, we use the 
quasiclassical Green's function theory.t!j~E!l The quasi- 
classical theory is a very powerful method, especiates 
for spatially wJiepiogeneous systems such as surfacesE^'EZl 
and vortices B3i23 Furthermore, one can easily treat a su- 
perconducting gap anisotropy as well as the Fermi surface 
anisotropy in the quasiclassical theory. We consider the 
transportlike Eilenberger equation for the quasiclassical 
Green's function 

?Aiuj r £)--iir( 3(w n ,T,Tk)_ i/(tw ni r,k) 
<?(^ n ,r,k)- «r^_. /t(i£( ^ rjk) _ 5 ( ia ,„, r , k ) 

(2.1) 

in a 2 x 2 matrix form (for even-parity superconductivity) , 
namely, 



iv F (k) ■ Vg(«w„,r,k) 

ioj n _ -A(r,k) 
A*(r,k) 



-ILJ„ 



$(iu; n ,r,k) =0. (2.2) 



The Eilenberger equation (2.2) is supplemented by the 
normalization condition 



-IT 2 !. 



(2.3) 



Here ui n = (2n + 1)ttT is the Matsubara frequency. The 
vector r = (x,y) is the center of mass coordinate, and 
the unit vector k represents the relative coordinate of the 
Cooper pair. The overbar denotes unit vectors. The com- 
mutator [A, B] = AB — BA. We assume the Fermi veloc- 
ity vp(k) is a function of k with reflecting the anisotropy 
of the Fermi surface. Since we consider an isolated single 
vortex in an extreme type II superconductor where the 
Ginzburg-Landau parameter k> 1, the vector potential 
can be neglected in Eq. (2.2). 

The Eilenberger equation in the matrix form (|]^) can 
be written down to the following equations, 

f w„ + r, 6) = A(r, %^„, r, 9), 



2 dr« 



v F (9) d 
2 dr« 



(2.4a) 



P(iu n ,T,6) = A*(r,%(tw n ,r,0), 

(2.4b) 
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v F (9)—g(ico n , r, 0) = A*(r, 0)/(iw n , r, 0) 
arii 

-A(r,e)/t( JWnjrj 0), 
which are supplemented by 

g(iuj n ,r,6) = [l-/(iw n ,r,fl)/t(iw„,r,e)] 1 / 2 , 
Re g(ioj n ,T,6) > 0. 

Here, k = (cos 0, sin 6*), 

v F (k) = (|v F (0)|cose(0),|v F (0)|sin6(0)) 
= (v F (0)cos9(0),u F (0) sin 6(0)), 



(2.4c) 



(2.5) 



(2.6) 



and the following coordinate system is taken: u = 
cos 6x + sin 0y , v = — sin Ox + cos 0y , and thus a point 
r = xx + yy is denoted as r = r\\u + r±v. The center of a 
vortex line is situated at the origin r = (0, 0). The angle 
0, i.e., the direction of k is measured from the a axis (or x 
axis) in the hexagonal plane of NbSc2 . If one considers a 
cylindrical Fermi surface with anisotropic Fermi velocity, 
then v F (k) = v F (0)k = (v F (0) cos 0, v F (0) sin 0) . 
The self-consistent equation is given by 



A(r, 0) = N 2irT 



UJ n >0 



d0' 
2^ 



p(0')V(9,0')f(iu Jn ,r,9'), 



(2.7) 



where Nq is the total density of states over the Fermi 
surface in the normal state. The 0-dependence of the 
density of states at the Fermi surface is represented by 



P(0) 



1 



JV |k- v F (k)| 



(2.8) 



which satisfies f (d9 / 2ir) p(9) — 1. We assume that the 
pairing interaction V(0,0') is separable, i.e., V(9,9') — 
vF(9)F(9'), where v is the strength of the pairing interac- 
tion and F(0) is a symmetry function, e.g., F{0) = cos 20 
for a d-wave pairing, F{9) = 1 for an isotropic s-wave 
pairing, etc. 

The pair potential is written as 



A(r,0) = A(r)F(0). 



(2.9) 



To o btain a sel f-consistent pair potential, we solve Eqs. 
( |2.4|) and (2/7) iteratively. This computation is per- 
formed after a method of Ref. |2^. In the calculation 
of the pap. potential, we adopt the so-called explosion 
methodfeaH to solve Eq. ([D|). 
The LDOS is evaluated from 



N(E,r)=N 



d0 
2^ 



p(0)Re g(iuj n — ► E + ir), r, 



^-p{9)N{E,v,0), 

Z7T 



(2.10) 



where rj (>0) is a small real constant. The value of 
r\ represents the effect of dilute impurities in a rough 



approximatior00 or other smearing effects!! To obtain 
g(iui n — > E+ir}, r, 0), we have to solve Eq. (|]J) for rj — iE 
instead of the Matsubara frequency u> n . While we suc- 
ceeded in this calculation in the vortex lattice case with 
the explosion method, a huge computer-running-timc for 
the numerical calculation was needed in this method.Ej 
In the case of the isolated single vortex, however, it is 
convenient to utilize ap-method of the Riccati equation 
developed by Schopohl.E§E§EJ The Riccati equation sim- 
plifies the numerical coroautation. 



The Riccati equationsc.il are given as 



v F (9) — a(w„,r,0) - A(r,( 
arii 



(2w„ + A*(r,0) a (a.„,r,0))aK,r,0) = 0, (2.11) 



v F (0)—b(io ni r,0) + A*(r, 

arii 



(2w„ + A(r, 0)b{u n , r, 0)) b(cu n , r, 0) = 0. (2.12) 



Equations ( |2.11 ) and (2.12) are obtained by substitut- 
ing the following parametrizationscJ into the Eilenberger 
equation (2.4), 



f = 



2a 



1 + ab' 



2b 



1 + ab' 



9 = 



(2.13) 



We solve Eqs. (2.11) and ( 2.12| ) independently along the 
r || -trajectory where r±_ is held constant. In the isolated 
single vorte x und er consideration, one can integrate Eqs. 
fl2.ll ) and ( 2.12 ) using solutions far from the vortex, 



a_ c 



V^ + |A(-oo,r x , 
A*(-oo,rj_, 



y/Lj* + |A(+oo,r±,( 
A(+oo,r ± ,0 



K>0) (2.14) 



as the initial values, respectively.c2l There, to remove the 
exploding solution, the integral about a is performed 
from ni = — oo, and about b from rn = -I-cxjMj We 
numerically integrate the first-order differential equa- 
tions (2.11) and ( p.l 1 ^ ) by the adaptive stepsize control 
Runge-Kutta method. O The Gre en's function g(iu> n — > 
E + ir ], r, 0) is ob tained from Eq. (2.12) if one solves Eqs. 
( |2.11 ) and ( 2.12| ) for rj — iE instead of u) n . When we 



solve Eqs. (2.11) and ( |2.12| ) for rj — iE, we use the self- 
consistently obtained pair potential A(r) which is calcu- 
lated beforehand. 

From now on, the density of states, energies, and 
lengths are measured in units of Nq, the uniform gap 
Ao at the temperature T = 0, and the coherence length 
Co = v F q/Aq (v fq = 1/N ), respectively. 
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III. PAIR POTENTIAL AND LOCAL DENSITY 
OF STATES 

Before going into technical details, we briefly explain 
our model and its parameter involved in connection with 
NbSe2, which is a typical type II s-wave superconductor. 
We assume the foll owing model of an anisotropic s-wave 
pairing in Eq. (2.S), 



F(6) = 1 + c A cos 60. 



(3.1) 



Here we again stress that the angle 8, i.e., the direction 
of k is measured from the a axis (or x axis) in the hexag- 
onal plane of NbSe2. Thus the parameter ca denotes 
the idewee of anisotropy in the superconducting energy 
gapc3~l£3 The case ca — corresponds to a conventional 
isotropic gap. 

The anisotropic s-wave gap is certainly suggested in 
NbSe2 from a scanningptunneling spectroscopy (STS) ex- 
periment at zero field.E3 The I-V tunneling spectrum, 
observed at the extreme low temperature T = 50 mK, 
indicates a substantial gap anisotropy (the gap ampli- 
tude with the averaged value 1.1 meV distributes from 
0.7 to 1.4 meV, see Fig. 1 in Ref. |23|), which is con- 
sistent with the density of states in the anisotropic s- 
wave gap case [Fig. |]J. It is seen from Fig. [I] that the 
gap edge distributes from E = 1 — ca to 1 + ca in the 
case ofrfche anisotropic gap. Then, the experimental data 
of STSE3 indicate that ca ~ 1/3. Similarly, a nuclear 
quadrupole resonance, NQR, experimented in NbSe2 sug- 
gests an anisotropic s-wave energy gap. The tempera- 
ture dependence of the spin-lattice relaxation rate 1 jT\ 
is well-fitted by an anisotropic energy gap model following 
HebeO with the value of a parameter S/A(0) ~ 1/3.C3 
Here the broadening in the gap edge, S/A(0), of Ref. Eq 
corresponds to S/tjq (0) ofRef.^. This parameter 6/ A (0) 
corresponds well to our parameter of the gap anisotropy, 
ca, because both parameters 6 and ca yield the broaden- 
ing in the gap edge. We set ca = 1/3 as a representative 
case in the following. 

In this paper, we restrict our attention to the gap 
anisotropy effect only, neglecting other effects, i.e., the 
vortex lattice effect and the effect of the anisotropic 
density of states at the Fermi surface, to clearly see 
how the energy gap anisotropy influences the LDOS. We 
calculate the LDOS in the isolated vortex case assum- 
ing an isotropic cylindrical Fermi surface (vp(k) || k, 
v F (9) =v F0 ). 



A. Pair potential 

In order to calculate the LDOS, we need the self- 
consistent pair potential obtained at the temperature, 
say, T = 0.1T C (T c is the superconducting transition tem- 
perature). The self-consistently obtained real-space vari- 
ation part of the pair potential, A(r), certainly exhibits 
a weak sixfold structure both in its phase and amplitude, 



which results from the anisotropic pairing, Eq. (| 
This behavior is similar to that of the d-wave case~Ea 
but sixfold symmetric here. In Fig. pi we show a contour 
plot of the amplitude of A(r). The amplitude |A(r)| is 
slightly suppressed in the x axis direction and its equiv- 
alent directions. As shown in Fig. ||, the spatial varia- 
tion of A(r) has weak anisotropy, but is almost circularly 
symmetric. However, the LDOS shows the characteristic 
sixfold symmetric structure as mentioned below. 



B. Local Density of States 

The LDOS calculated using the self-consistent pair po- 
tential has almost the same structure, except for the 
length scale for its spread, as that calculated using 
a test-potential A(r) = A(T) tanh(r/£) exp(i</>) does, 
where A(T) is the uniform gap at the temperature T, 
£ = t>Fo/A(T), and the cylindrical coordinate system 
r = (r, <fi) is taken. That is, the LDOS does not so af- 
fected by the weak sixfold symmetric spatial structure of 
the real-space variation part of the pair potential, A(&L 
We have seen the same situation also in the d-wave case.E2l 
It means that a calculated sixfold or fourfold structure 
of the LDOS directly results from the fc-space variation 
part of the pair potential, F (9). 

In Fig. ||, we show the LDOS N(E, r) for several en- 
ergies E in the case ca = 1/3, calculated by using the 
self-consistently obtained pair potential. It is seen from 
Fig. ||(a) that the sixfold star centered at the vortex cen- 
ter is oriented away from the x axis by 30° for E = 0. 
Next it is seen from Fig. ||(b) that at the intermediate 
energy each ray splits into two parallel rays, keeping its 
direction. This characteristic featuce—was precisely ob- 
served in the experiment by HessOE3 With increasing 
the energy E further, the sixfold star becomes a more 
extended one, and its orientation rotates by 30° as seen 
from Fig. ||(c). Note that in Fig. ||(c) the head of each 
ray splits in two. It coincides with an experimental re- 
sult (see the STM image for 0.48 mV in Fig. 1 of Ref. 
^l| ). In this way, the anisotropic s-wave gap model well 
reproduces the experimental features mentioned in Sec. 
[j| (1) the sixfold star shape, (2) the 30° rotation, and 
especially (3) the split parallel ray structure at the inter- 
mediate energy. We refer to Fig. 1 in Ref. |3l] where the 
density plots of the LDOS compared with the experimen- 
tal data are displayed, which is complimentary to Fig. || 
in the present paper. 

Another way to examine the quasiparticle excitations 
in the vortex states is to see how the spectrum evolves 
along radial lines from the vortex center. We show the 
spectral evolutions along the radial lines for 30° in Fig. 
|(a), 15° in |(b), and 0° in |(c) from the x axis. The 
zero-bias peak splits into several peaks in each spectral 
evolution. Cross sections of each spectral evolution at 
r = 1 (r = \J x 2 + y 2 ) are shown in Fig. [5] to provide the 
identification of each ridge in Fig. [|. 
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In the calculation of Figs. || and f|, the smearing fac- 
tor is chosen as r\ — 0.03, which well reproduces the 
STM experimental data. It corresponds to the solid lines 
of Fig. where the peaks are labeled a-e. The case 
with smaller smearing effect (rj — 0.001) is represented 
by the dashed lines in Fig. where the spectrum has 
the sharp peaks labeled as A-E. (The structure of these 
peaks is discussed in the next section.) As shown in Fig. 
||, by increasing the smearing effect, the spectrum of the 
dashed line (77 = 0.001) is reduced to that of the solid 
line (77 = 0.03), and reproduces the STM experimental 
data. It seems that the LDOS actually observed in STM 
experiments is somewhat,-smeared due to impuritiesuj'tS 
or other smearing effects.c3 

In Fig. ^(a) (the 30° direction), there exist one peak 
at E = and three pairs of peaks. The peak at E = 
in Fig. [|(a) [the e peak in Fig. |](a)] corresponds to the 
ray which extends in the 30° direction in Fig. 0(a). This 
peak is referred to as the inner peak in Refs. ^ and |2^. 
This inner peak [the e peak] corresponds to also the split 
parallel ray in Fig. |5|(b ) and the head of the ray which 
splits into two in FigTWc). The inner e peak is, therefore, 
sensitive to the angle of the radial line, and splits in a 
pair of peaks with the variation of the angle [see Figs, 
^(b) and|^(c)]. On the other hand, the most inside pair 
of peaks in Fig. |](a) [the 5 peak in Figs. |J and H is not 
sensitive to_the angle. This peak is referred to as the 
outer peakocj The behavior of the calculated inner and 
outer peaks well coincide with the experimental result 
(the experimental finding (4) in Sec. |j). The positions 
of the outer 5 and inner e peaks as a function of r are 
compared with the experimental data in Fig. 3 of Ref. 

m 

Outside the inner e and outer 8 peaks, extra peaks ap- 
pear in each calculated spectral evolution [the a, /3, and 
7 peaks in Figs. || and The result of the calculation 
shows that the extra peaks are relatively sensitive to the 
angle of the radial line. The existence of the extra peaks 
is characteristic of the gap anisotropy effect. The peaks 
a and f3 merge into the upper edge of the energy gap, 
1 + ca, far from the vortex. These extra peaks have not 
been noted in experimental data so far. While each peak 
cannot be clearly identified in experimental data yet, it 
seems that there is_at least one new line outside the outer 
peak in the dataJ£j It is expected for future experiments 
to definitely identify the extra peaks. 

The dependence of the LDOS on the angle of the radial 
line is important, because it gives a detailed information 
on the gap anisotropy. To see it, we show in Fig. ^ a 
spectral evolution from the angle 0° to 30° along a circle 
whose radius r = 1. From this, we can see how each peak 
moves, and joins up the others with the variation of the 
angle. As mentioned above, the e peak (that is, inner 
peak) is sensitive to the angle cj> of the radial line, and 
the 6 peak (outer peak) is insensitive to (j>. The e peak is 
located at E = for = 30°. When <p deviates from 30°, 
the peak splits into two which are positive and negative 
energy peaks. With decreasing <p to 0°, the energy im- 



position of the £ peak increases. As for the peaks a, /3, 
and 7, with decreasing from 30° to 0°, the iS-position 
decreases for the 7 peak, increases for the (3 peak, and 
is insensitive for the a peak. The peaks (3 and 7 overlap 
each other for tj> = 30°, and the peaks a and (3 overlap 
each other for = 0° (see also Fig. ||). Here, we should 
mention the behavior of the 7 peak at <fi ~ 0°. In Fig. ||, 
the 7 peak seems to join up the angle-insensitive 6 peak 
near 0°, that is, the 7 peak is buried in the (5 peak in Figs. 
||(c) and ||(c) (the 0° direction) . Such a behavior of the 
7 peak intimately relates to the value of the anisotropic 
gap parameter, ca- The above behavior of 7 is that of 
the case ca = 1/3. According as ca increases further, the 
position of the 7 peak at 0° shifts to the higher energy 
side [see a spectral evolution shown in Fig. 3(a) of Ref. |3l] 
(the 0° direction), where ca is set to 1/2 and we can see 
that a peak line, which corresponds to the present peak 7 
(not denoted explicitly in that figure), evolves away from 
the 5 peak line]. 



IV. QUASIPARTICLE TRAJECTORIES 

In this section, we interpret the behavior of the quasi- 
particle bounded around a vortex in terms of the quasi- 
classical picture. 



A. Direction-dependent Local Density of States 

In the quasiclassical approximation, the equations are 
independently given for each direction of k. The Eilen- 
berger equation (or the Riccati equation) for a direction 
k is independent of those for the other directions. The 
direction-depende nt lo cal density of states N(E, r, 9) in- 
troduced in Eq. (2.10) is obtained from the solution of 
the equation for the direction k = (cos 9, sin 9). The 
LDOS N(E, r) is calculated by integrating the direction- 
dependent LDOS N(E, r, 9) over 9. In an isolated vortex 
state, the structui^-of-JV (E, r, 9) was previously investi- 
gated analyticallyQ'EjO and numerically!!^ According to 
the results of these investigations, N(E, r, 9) has the fol- 
lowing structure .for low energies below Ao in the isolated 
single vortexlalijO (Here, remind ourselves of the nota- 
tion: r = ix + yy — riiii + r±y] u = cosf?x + sin#y, 
v = — sinflx + cos6*y.) (i) N(E,r,9) as a function of 
r = (ri|,ri) vanishes everywhere except on a straight 
line along which r± — const. = rj_(E). This straight 
line and r±(E) are referred to as "quasiparticle path" 
and "impact parameter," respectively, (ii) Along the 
line r_L = r±(E), N(E,r,9) has a single maximum at 
rn = and decreases exponentially for r\\ — ► ±00. (iii) 
The impact parameter r± (E) is a monotonically increas- 
ing function of E. One defines E(r±) as the energy level 
of the state on the quasiparticle path with the impact pa- 
rameter rj_. In extreme type II superconductors where 
/t>l, E(r±) is determined by the minimum value of the 
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amplitude of the pair potential on the quasiparticle path 
r_i_ = r±(E). For the low energy levels, E(r±) is given by 
E(r±)=sgn(r±)\A(r\\ = 0,r±)\ in a good approximation. 

On the basis of the above properties (i)-(iii) of the 
direction-djapendenLLDOS N(E, r, 0) f studied by Kramer 
and Peschja Klein Jl3 and Ullah et aZ.jLZl we interpret our 
result of the preceding section as follows. 

For simplicity, we concentrate our attention to Eq. 
(2.11) as a representative. Dividing Eq. (2.11) by F(6), 
we rewrite this equation as 



1 



F(6) dr 



-a(w„,r,0) - A(r) 



+ + A*(r)a( W „, r, 9)) a{u n , r, 6) = 0. (4.1) 

In the case of the isotropic s-wave pairing (F(9) = 1), 
N(E, r, 9) at a fixed energy has the identical structure for 
each direction 9 [the items (i), (ii), and (hi)]. Then the 
LDOS N(E,r), obtained by integrating N(E,r,0) over 
9, exhibits a "ring" shaped structured in the real space. 
The impact parameter is the radius of the ring. 

In the case of an anisotropic pairing, the situation is 
chan ged because of the te rms which include F(8) in Eq. 
(4.1). According to Eq. (4.1), both the length scale in 
the r || -direction and the energy scale vary with 9, but 
otherwise the form of the equation is same as that of the 
isotropic s-wave case. For the direction 9 where F(9) 
is suppressed, the length of the spreading of N(E, r, 9) 
along the quasiparticle path [note the items (i) and (ii)] 
becomes large. For the same 9, the effective energy be- 
comes large and then the impact parameter becomes far 
from the vortex center [note the item (hi)] . 



B. Interpretation on the LDOS around a vortex 

We show the partly integrated N(E,r,ff) in Fig. |t], 
where the integration is done from 9 = —30° to 30°, 
and its schematic figure in Fig. |^, for the pairing of Eq. 
(3.1) where ca = 1/3. Here, to clarify the structure of 
the LDOS, a small smearing parameter (rj = 0.001) is 
adopted. The peak lines shown in Fig. ^ are composed 
of the quasiparticle paths of each direction 9 described 
above. These peak lines can be interpreted as the flows 
of quasiparticles shown in Fig. ||. It is noted that the 
trajectories 1 and 2 appear, because F(9) is finite at 
9 = —30° and 30°, i.e., the impact parameter is finite 
at these angles. If F(9) has a node, i.e., ca = 1, the 
impact parameter is infinitely far from the vortexj-eenter 
for the quasiparticle path of the node direction,^ and 
the trajectories 1 and 2 disappear. In the bound states, 
the quasiparticles flow along these trajectories. We call 
it "quasiparticle trajectory." The whole state at a fixed 
energy is composed of such flows of quasiparticles along 
the quasiparticle trajectories, while the individual quasi- 
particle paths of each direction 9 [the items (i) and (ii)] 
could be considered to be the Andreev reflections. 



We show in Fig. | the LDOS N(E, r) obtained by in- 
tegrating the direction-dependent LDOS N(E, r, 9) over 
all 9. A schematic figure which corresponds to Fig. ^| is 
shown in Fig. [h]. The peaks which the radial lines cross 
are labeled A-E there. When the energy E elevates, the 
scale of the trajectory in Fig. 10 increases with keeping 
its structure fixed. Therefore, the trajectory has one-to- 
one correspondence to the peak of the spectrum of Figs. 

The peaks A-E of Fig. [lO] precisely correspond to 
those of Fig. |j| These peaks are smeared to appear as a- 
e peaks in Fig. | (and thus in Figs. | and |). The LDOS 
actually observed in STM experiments is not that shown 
in Fig. ^ itself, but somejwhat smeared one [Figs. ^, 1 
and 0] due to impurities! 1 6 H 18 1 or other smearing effects.EHI 
Roughly speaking, the peaks A, B, C, Dl, and E corre- 
spond to the peaks a, (3, 7, <5, and e, respectively. 

The trajectory of Fig. |lfj helps us to facilitate an un- 
derstanding of the rich structure of the LDOS. The tra- 
jectories B and C cross each other at the angle 4> = 30° 
from the x-axis in Fig. ful Then, the peaks B and C (i.e., 
/3 and 7) overlap each other in Figs. |(a) and |](a). The 
cross of the trajectories A and B at 4> = 0° in Fig. |To| 
corresponds to the overlap of the peaks A and B (i.e., a 
and /3) in Figs. ||(c) and ||(c). When <fi varies from 30° 
to 0°, the trajectories C and Dl cross each other in Fig. 
|Io| , where ca = 1/3. It corresponds to the result that 
the peaks 7 and 5 interchange their positions between 
Figs. ||(b) and ||(c). However, this behavior of 7 and S 
depends on the anisotrop_ic gap parameter ca as men- 
tioned at the end of Sec. 
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In the case of large ca , the 
trajectories C and Dl does not cross for 0° < < 30° 
in Fig. [l(]. Even at <f> — 0° , the trajectory Dl is located 
farther from the vortex center than the trajectory C, for 
large ca- Then, the peak C (i.e., 7) is located at higher 
energy than the peak Dl (i.e., 8) in the spectrum of Figs. 
||(c) and ||, for large ca- As seen in Fig. |[ the peak D2 
tends to be buried in the other peaks, due to the smear- 
ing effects. However, if the experiment is performed for 
the weak smearing case, the peak D2 should be observed 
as a small peak, which splits from the peak Dl (i.e., S) at 
4> = 30° and approaches the peak E (i.e., e) with decreas- 
ing tf> to 0°. This D2 peak seems to be easily observed 
for the angle 0° < <fi < 10°. We detect a small indication 
of the D2 peak for this angle region, if Fig. || is enlarged 
at <f> ~ 0°. The trajectories Dl and E (i.e., 5 and e) cor- 
responds to the trajectories 1 and 2 of Fig. ^, which is 
related to the lower edge of the anisotropic energy gap. 
Therefore, these trajectories disappear for the higher en- 
ergy, (1 - ca) < E < (1 + ca). The peaks Dl and E (i.e., 
6 and e) merge into the lower edge of the energy gap at 
E = 1 — ca far from the vortex. 



C. Flows of quasiparticles around a vortex 

The flows of the quasiparticles mentioned above are 
quantitatively represented by the following quantity, 
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1(25, r) = (* ^ p (6)v F (0)N(E, r, 9), (4.2) 

which we tentatively call "directional local density of 
states." This directional LDOS corresponds to a quan- 
tity obtained by integrating-, "spectral current density" 
introduced by Rainer et al£$ over 9 (or p/ in Ref. |39[ ). 
The total current density arouiad a vortex is composed 
of the spectral current density.ta In Figs. |ll](a), |Tl|(b), 
|ll](c), and |l|(d), we show the directional LDOS 1(25, r) 
calculated for E = 0.2, 1.2, 1.4, and 1.6, respectively. 
Here 1(E, r) is calculated under the condition considered 
in this section, i.e., under the anisotropic gap and the 
isotropic cylindrical Fermi surface. It is seen from Fig. 
|ll|(a) that the flow of the quasiparticle exhibits a sixfold 
anisotropy resulting from the sixfold LDOS of the bound 
states (Fig. ||(b) and thus Fig. [H]). Now, it is of interest 
to note the flow with an energy near the upper gap edge, 
E = l + c A 1.3). Comparing Figs. |T|( b ) and 0(d), we 
can see that the quasiparticles above and below the up- 
per gap edge flow each other in reverse directions except 
in the vicinity of the vortex center. It certainly coincides 
with a result of an analysis, based on the Bogoliubov- 
de Gennes (BdG) equation.liJ This feature should not be 
influenced by the gap anisotropy. 

V. SUMMARY AND DISCUSSIONS 

The LDOS around an isolated single vortex is studied 
within the framework of the quasiclassical theory. We 
consider the effect of the anisotropy of the superconduct- 
ing energy g ap. Assuming the anisotropic s-wave energy 
gap in Eq. ( |3.l| ), we succeed in theoretically reproduc- 
ing the characteristic structure of the LDOS observed in 
STM experiments; the observed features, i.e., the items 
(l)-(4) for NbSe2 listed in Sec. [j| are well described in 
terms of the anisotropic gap model. We point out the ex- 
istence of the missing peaks (a, (3, and 7) at the higher 
energy side in the spectral evolution shown in Figs. 0-pl 
which is expected to be looked for in a future experiment. 
We also notice the further splitting of the observed broad 
peaks as shown, for example, in Fig. 5(b) (<5 — > Dl and 
D2). These predictions, which reflect the gap anisotropy, 
may be checked by using a purer sample at lower temper- 
atures, because smearing effects, due to lattice defects or 
thermal broadening, mask the fine details. We attempt to 
interpret the calculated LDOS in terms of the quasiparti- 
cle trajectory. This enables us to thoroughly understand 
the STM results and the internal electronic structure of 
the vortex. In this paper, the value of our parameter is 
chosen appropriate for NbSe 2 . However, the essence of 
the obtained results should be applicable to other type II 
superconductors in general although the degree of the gap 
anisotropy ca and the symmetry of F(9) will be different 
in each case. Even in other anisotropic superconductors, 
the explanation in terms of the quasiparticle trajectory 



would be helpful to an understanding of the internal elec- 
tronic structure of vortices. 



A. Comparison with other theories and effects of the 
vortex lattice 

Let us comment on prior works which are connected 
with the star-shaped LDOS observed in NbSe2- On, the 
basis of a sixfold perturbation, Gygi and SchluterEll ex- 
plained that the lower and higher energy stars observed 
by STM were interpreted as bonding or antibonding 
states. The STM results (1) and (2) listed in Sec. | were 
able to be explained by this perturbation scheme. They 
adopted a sixfold crystal lattice potential in NbSe2 as, 
the perturbation. Recently, Zhu, Zhang, and SigristE3 
investigated the effect of the underlying crystal lattice 
by means of a non-perturbation method, i.e., a method 
of diagonalizing a tight-binding BdG Hamiltonian in a 
discrete square lattice, where the crystal lattice poten- 
tial, i.e., the band structure is determined a priori. This 
method supplements the perturbation theory of Ref. |2^: 
the absolute orientation of the star relative to the under- 
lying crystal lattice was determined.E3 

By this non-perturbation approach, also a gradual ro- 
tation of the star-shaped LDOS was obtained in the in- 
termediate energy region.E3 Nevertheless, it is not yet 
clear whether the crystal lattice effect is able to repro- 
duce the remaining experimental findings (3) and (4), 
i.e., the split parallel ray structure and the behavior of 
peaks in the spectral evolutions. The model used in Ref. 

is the discrete lattice model, and therefore it is impos- 
sible to obtain detailed spectra, e.g., spectral evolutions 
along radial lines, due to the discreteness. Hence, it is 
desired to treat the crystal lattice potential effect with a 
non-perturbation method in the continuum limit. 

Now, the crystal lattice potential determines the band 
structure, and influences the structure of the Fermi sur- 
face. The effect of the crystal lattice potential should ap- 
pear as the anisotropy of the Fermi surface. In our frame- 
work, the anisotropy of the Fermi surface is taken into ac- 
count by assuming an anisotropic density of states at the 
Ferm i surface, p{9), which appears in the 9- integral of Eq. 
( |2.10| ), and the anisotropic Fermi velocity vp(9), which 
appears in the Eilenberger (or Riccati) equations. The 
experimental findings (l)-(4) can be reproduced quali- 
tatively, if we introduce a large anisotropy in v-p(9). Its 
details and the comparison of them to the present results, 
from the gap anisotropy effect are discussed elsewhere.ES 
Gygi and Schliiter considered also the effect of nearest- 
neighbor vortices, i.e, that of the vortex lattice. El They 
adopted a sixfold anisotropy of the vector potential as the 
vortex lattice effect, and treated it as the perturbation. 
However, the periodicity of the pair potential is also im- 
portant as the effect of the vortex latticeEElo In extreme 
type II superconductors such as NbSe2 where k 3> 1, 
the periodicity of the pair potential is expected to have 
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stronger effects upon the structure of the LDOS than the 
anisotropy of the vector potential does. We find in Ref. 

that the effect of the periodicity gives a characteristic 
sixfold structure to the LDOS. 

This structure of the LDOS which results from the pe- 
riodicity of the pair potential appears only at high mag- 
netic fields such as 1 T for the material parameters ap- 
propriate to NbSe2, where the vortex core regions sub- 
stantially overlap each other l3 At a lower magnetic field 
such as 0.1 T, the calculated LDOS reduces to the al- 
most circular structure. On the other hand, the LDOS 
observed in a STM experiment exhibits the star-shaped 
structure in spite of a low field 0.025 T (see Fig. 12 in 
Ref. [53|). Therefore, in the case of NbSe2 at low mag- 
netic fields, we need to consider the effects of anisotropy 
other than the vortex lattice effect in order to explain the 
star-shaped LDOS. Both the vortex lattice effect and the 
anisotropic superconducting gap one are important for 
the star-shaped LDOS observed in NbSe 2 at high mag- 
netic fields. We expect a future STM experiment to be 
performed on isotropic superconducting compounds or 
metals to clarify the vortex lattice effect and confirm pre- 
dictions of Ref. [3l. 

In STM experiments on NbSe2, one of the directions 
of nearest-neighbor vortices coiad das— .with .the a axis 
(see the literature by Hess et alMaBlljMki or Ren- 
ner et al.B), except for extreme low fields.E3 This exper- 
imental fact gives evidence of a correlation of the vor- 
tex lattice with the underlying crystal lattice of NbSe2. 
It was recently found that in d-wave superconductors, 
higher-order (nonlocal correction) terms in the Ginzburg- 
Landau equation, which reflect the fourfold symmetric 
property of the <i-wave pairing .-give rise to a preferred 
direction of the vortex lattice.E3l£3 In NbSe2, the sixfold 
anisotropic pairing, Eq. fl3.ip , is expected to give rise to 
the same correlation as the d-w&ve pairing does, and it 
may be the origin of the experimental fact mentioned 
above. A possibility of the correlation of the vortex lat- 
tice with the underlying crystal lattice was recently re- 
ported also in a high T c cuprate. l3'E3 

B. Beyond the quasiclassical approach 

We mention the LDOS around a vortex in high T c 
cuprates. It seems from various experiments that high 
T c material is a d-wave superconductorEJ A fourfold 
structure of the LDOS is predicted in d-wave supercon- 
ductors byJthcaictical studies based on the quasiclassi- 
cal theoryoESEEl The origin of this fourfold structure 
is same as that discussed in the present paper for the 
gap anisotropy. Recently, Maggio-Aprile et al. observed 
tunneling spectra around yorliEes in a high T c cuprate, 
YBa2Cu 3 07_,5, with STmMBH However, the spectro- 
scopic images of STM have not exhibited any sign of a 
fourfold structure yet. We expect further detailed experi- 
ments to observe the fourfold symmetric LDOS structure. 



When we consider the high T c materials, the quantum 
effects should be taken into account. The quasiclassi- 
cal theory is certainly valid only in systems where the 
atomic scale spatial variation of the Green's function can 
be neglected with respect to the coherence length scale 
one.ESH The effects neglected in the quasiclassical the- 
ory can be important in the case of the high T c cuprate; 
the quantization of energy levels of the bound states can- 
not be treated by the quasiclassical theory, and while it 
is possible in the quasiclassical approximation to divide 
the equation into individual equations for each direction 
of k, it is impossible in the quantum-mechanical limit. 
Although we expect the fourfold structure of the LDOS 
should be observed in future experiments, the above ef- 
fects may change the situation in the case of the high 
T c cuprate. It is certainly desired on the theoretical side 
that a fully quantum-mechanical approach clears up this 
problem in future. 

As for the fully quantum-mechanical approach, it is 
needed to solve the BdG equation without quasiclassical 
approximations. The BdG equation cannot be written 
in a local form in the case of an anisotropic pairing, and 
therefore it is difficult to treat this equation in the con- 
tinuum limit. One of the possible approaches to this 
problem is the method of diagonal i zi ri g_ a BdG Hamil- 
tonian for a specific lattice model. In the lattice 
model, however, the atomic scale variation of wave func- 
tions among the lattice points is uncertain. 

In most superconductors (£ S> 1/&f), the atomic scale 
variation of the wave function is a redundant informa- 
tion and can usually be neglected on the basis of the 
quasiclassical theory. On the other hand, in the high T c 
cuprate superconductors, fcp£ — £f/Ao ~ 1 (the Fermi 
wave-number and energy are fcp and £f ; respectively) ,l3 
and therefore the atomic scale variation and the quanti- 
zation of bound states in a vortex may be crucial for the 
electronic structure around the vortex in the cuprates. 

The high T c cuprate is certainly the only supercon- 
ductor possessed of a possibility of an experimentally de- 
tectable quantization in the vortex bound states. Accord- 
ing to Ref. [j], a substantial energy quantization (of the 
order of Ap/ep ~ 10 K) is expected to exist in the high 
T c cuprate. However, to the present author's knowledge, 
the system considered in Ref. |l| is an isotropic s-wave 
superconductor and the mechanism of the quantization 
in the case of anisotropic pairing is not yet understood. 
In case of gap node due to anisotropic pairing, it is ex- 
pected that the separation of the energy levels becomes 
small. Further experiments, which, e.g., investigate spa- 
tial variation of this quantized bound states in the high 
T c cuprates with STM and [then, compare its result with 
the quasiclassical predictioroHH in order to clarify how 
the quantum effects mentioned above modify the vor- 
tex bound states, are the need for alternative theoretical 
studies of the vortex bound states. 
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C. Concluding remarks 

The electronic structure of vortices in a compojAiwi, 
LuNi 2 B 2 C, was quite recently investigated by STM.E3'E2l 
Although no conductance peaks related to localized 
quasiparticle states in the vortex core are observed in the 
experiment, due to a short mean free path (of the order of 
the coherence length), and thermal broadening effects at 
4.2 K (T c « 16 K)jB a rich (maybe fourfold) structure 
of the LDOS such as that discussed in the present pa- 
per is expected to be detected in STM spectra by lower- 
ing the temperature and decreasing impurities or defects. 
If an anisotropic bound states around a single vortex is 
observed, it should suggest an anisotropy of the pair- 
ing in this compound. The direction (in the A-space) in 
which the superconducting gap is suppressed corresponds 
to that (in the real space) of a ray of the LDOS at zero 
bias. 

Finally, low-temperature STM is the unique experi- 
mental method which has the ability not only to image 
the distribution of the vortex lattice, but also to probe 
the electronic structure of individual vortices. We ex- 
pect future STM experiments to be performed in vortex 
states on various superconductors such as organic con- 
ductors, high T c cuprates, heavy Fermion superconduc- 
tors (e.g., UPt 3 ), and a recently discovered Jipn-copper- 
layered perovskite superconductot^-S^RuOj^ which has 
nearly cylindrical Fermi surfacealjQ and a possibility 
that lEiLjadd-parity superconductivity would be realized 
in itO'La The information on the vortex bound states 
available from STM spectra can be one of clues to the 
pairing. The low-tcmperature STM experiments deserve 
a great deal of attention. 
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FIG. 11. The directional LDOS I(E, r) (77 = 0.03) for the 
energies E = 0.2 (a), 1.2 (b), 1.4 (c), and 1.6 (d). The arrows 
in the figures represent only the directions of I(E,r). 



FIG. 2. Contour plot of the amplitude of the real-space 
variation part of the pair potential, |A(r)|. From the center, 
0.1, 0.2, 0.9. The temperature is T = 0.1T C , and the 
degree of the gap anisotropy is ca = 1/3. 



FIG. 3. The LDOS N(E, r) (77 = 0.03) calculated for the 
energies E = (a), 0.2 (b), and 0.4 (c). Large peaks in the 
vicinity of the vortex center are truncated in the figures (a) 
and (b). 

FIG. 4. Spectral evolutions N(E, r) (77 = 0.03) along ra- 
dial lines for 30° (a), 15° (b), and 0° (c) from the x axis. The 
zero-bias peak is truncated in the figures. The peak lines in 
the spectra are labeled a-e. 

FIG. 5. Cross sections of the spectral evolutions (Fig. 4) 
at the distance from the vortex center, r = 1. The directions 
of each radial line are 30° (a), 15° (b), and 0° (c) from the 
x axis. The peaks in the spectra are labeled A-E for the 
dashed line spectra (77 = 0.001) and a-e for the solid line 
spectra (77 = 0.03). The labels a-e correspond to those of 
Fig. 4. 

FIG. 6. Spectral evolution N(E, r) (77 = 0.03) from the 
angle 0° to 30° along a circle whose radius r = 1. The cen- 
ter of this circle is situated at the vertex center. The peaks 
labeled as a-e correspond to those of Figs. 4 and 5. 

FIG. 7. The direction-dependent LDOS N(E, r, 9) partly 
integrated from 9 = -30° to 30°, where E = 0.5, 77 = 0.001, 
and 6£o x 6£o is shown in the real space. 

FIG. 8. Schematic flow trajectories of quasiparticles with 
an energy < E < (1 — ca). These trajectories correspond 
to those shown in Fig. 7. When (1 — ca) < E < (1 + ca), 
the trajectories 1 and 2 disappear and only the trajectory 3 
is alive. 

FIG. 9. The LDOS N(E, r) which is obtained by integrat- 
ing the direction-dependent LDOS N(E, r, 9) over 9, where 
E = 0.5, 77 = 0.001, and 6£o x 6£o is shown. 

FIG. 10. Schematic figure of the LDOS N(E, r) for an 
energy < E < (1 — ca). Points A-E correspond to the 
peaks of the dashed line spectra in Fig. 5. 



FIG. 1. The density of states N(E) at zero field, where 
smearing parameter 77 = 0.03 (solid line) or 77 = 0.001 (dashed 
line). It is calculated for the isotropic cylindrical Fermi sur- 
face, and the degree of the gap anisotropy is ca = 1/3. 
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